#include <cmath>
#include <cassert>
#include <stdlib.h>
#include "Complex.h"
#include "SectionIntegrals.h"

/*
 * C version of r1mach and d1mach from core (netlib)
 * specialized for IEEE arithmetic
 * 
 * MACHINE CONSTANTS (s for single, d for double)
 * {S|D}1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
 * {S|D}1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
 * {S|D}1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
 * {S|D}1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
 * {S|D}1MACH(5) = LOG10(B)
 */
// #define M_PI 3.14159265358979323846264338328
// #define M_LN2 0.693147180559945309417232121458
static double	mach[5] = {2.2250738585072014e-308,
			   1.7976931348623157e+308,
			   1.1102230246251565e-16,
			   2.2204460492503131e-16,
			   3.0102999566398120e-01};

static const double s11r[21] = {
  0.0448875760891932036595562553276,0.0278480909574822965157922173757,
  0.00394490790249120295818107628687,-0.00157697939158619172562804651751,
  -0.0000886578217796691901712579357311,0.0000301708056772263120428135787035,
  9.521839632337438230089618156e-7,-3.00028307455805582080773625835e-7,
  -6.14917009583473496433650831019e-9,1.85133588988085286010092653662e-9,
  2.67848449041765751590373973224e-11,-7.82394575359355297437491915705e-12,
  -8.44240072511090922609176843848e-14,2.41333276776166240844516922196e-14,
  2.02015531985181413114834031833e-16,-5.68171271075270422851146478874e-17,
  -3.80082421064644521052871349836e-19,1.05551739229841670238163200361e-19,
  5.7758422925275435667221605993e-22,-1.58774695838716531303310462626e-22,
  -7.24181766014636685673730787292e-25
};

static const double s11i[21] = {
  0.100116671557942715638078149123,0.0429600096728215971268270800599,
  -0.00799014859477407505770275088389,-0.000664114111384495427035329182866,
  0.000240714510952202000864758517061,9.89085259369337382687437812294e-6,
  -3.22040860178194578481012477174e-6,-8.08401148192350365282200249069e-8,
  2.48351290049260966544658921605e-8,4.24154988067028660399867468349e-10,
  -1.25611378629704490237955971836e-10,-1.56053077919196502557674988724e-12,
  4.50565044006801278137904597946e-13,4.2641179237225098728291226479e-15,
  -1.2084245714879456268965803807e-15,-9.01338537885038989528688031325e-18,
  2.5180796700698002962991581923e-18,1.51955263898294940481729370636e-20,
  -4.19737873024216866691628952458e-21,-2.092488792285595339755624521e-23,
  5.72708467031136321701747126611e-24
};

static const double s12r[21] = {
  -0.376145877558191778393359413441,0.0775244431850198578126067647425,
  0.0120396593748540634695397747695,-0.00385683684390247509721340352427,
  -0.000232359275790231209370627606991,0.0000697318379146209092637310696007,
  2.32354473986257272021507575389e-6,-6.71692140309360615694979580992e-7,
  -1.43946361256617673523038166877e-8,4.06087820907414336567714443732e-9,
  6.10183339004616075548375321861e-11,-1.69196418769523832825063863136e-11,
  -1.88669746820541798989965091628e-13,5.16473095452962111184823547686e-14,
  4.45066881692009291504139737861e-16,-1.20625107617859803735741992452e-16,
  -8.28193837331508300767103116139e-19,2.22680015825230528892642524445e-19,
  1.24755889505424049389100515561e-21,-3.33254971913153176741833960484e-22,
  -1.55307002839777371508497520751e-24
};

static const double s12i[21] = {
  0.0527472790869782317601048210983,0.00823962722148093961886198320927,
  -0.0205185842051817330153151013327,-0.00184683218270819613487368071941,
  0.000569681886932212757533488372406,0.0000248774530818801164177266528608,
  -7.31121019876580624171992432347e-6,-1.92744564223806538367454388776e-7,
  5.49794278719049727550379096876e-8,9.78237385539447442446850072421e-10,
  -2.7341624177723508216430132999e-10,-3.51839815887772323640101921381e-12,
  9.68934411607055794052256859665e-13,9.45703963505047353201918875825e-15,
  -2.57516976113400217760868402425e-15,-1.97419921753098238455550504742e-17,
  5.32820017906655555903355375475e-18,3.29581793797656865402793252539e-20,
  -8.83137325823594007269279476114e-21,-4.50279718100548728336329365981e-23,
  1.19941679774924468309434420379e-23
};

static const double m12r[21] = {
  0.148523151773238914750879360089,-0.0117856118001224048185631301904,
  -0.00248887208039014371691400683052,0.000250045060357076469386198883676,
  0.0000227217776065076434637230864113,-2.48764935230787745662127026799e-6,
  -1.32138506847814502856384193414e-7,1.50966754393693942843767293542e-8,
  5.3472999553162661403204445045e-10,-6.26136041009708550772228055719e-11,
  -1.59574066624737000616598104732e-12,1.89788785691219687197167013023e-13,
  3.66030609080549274006207730375e-15,-4.39955659500182569051978906011e-16,
  -6.65848768159000092224193226014e-18,8.06343127453005031535923212263e-19,
  9.84397490339224661524630997726e-21,-1.19869887155210161836484730378e-21,
  -1.20634550494837590549640883469e-23,1.47512193662595435067359954287e-24,
  1.24549093756962710863096766634e-26
};

static const double m12i[21] = {
  -0.0454399665519585306943416687117,-0.0210517666740874019203591488894,
  0.00194647501081621201871675259482,0.000253466068123907163346571754613,
  -0.0000268083453427538717591876419304,-1.82138740336918117478832696004e-6,
  2.04357511048425337951376869602e-7,8.75944656915074206478854298947e-9,
  -1.01466837126303146739791005703e-9,-3.02573132377805421636557302451e-11,
  3.57358222114420372764650037191e-12,7.88121312149152771558608913996e-14,
  -9.42758576193708862552405242331e-15,-1.60439904050827900099939709069e-16,
  1.93624791035947590366500765061e-17,2.62394448214143482490534256935e-19,
  -3.18700789496399461681365308408e-20,-3.52400207248027768109209530864e-22,
  4.30074555255053206057921088056e-23,3.95655079023456015736315286131e-25,
  -4.84642137915095135859812028886e-26
};

/*
 * evaluate a chebyshev series
 * adapted from fortran csevl
 */
static double
csevl( const double x, const double *cs, unsigned int n )
{
  double b2 = 0, b1 = 0, b0 = 0;
  const double twox = 2 * x;

  while( n-- ){
    b2 = b1;
    b1 = b0;
    b0 = twox * b1 - b2 + cs[n];
  }

  return ( b0 - b2 )/2;
}

/*
 * from the original fortran inits
 * april 1977 version.  w. fullerton, c3, los alamos scientific lab.
 *
 * initialize the orthogonal series so that inits is the number of terms
 * needed to insure the error is no larger than eta.  ordinarily, eta
 * will be chosen to be one-tenth machine precision.
 */
static unsigned int
inits( const double *series, unsigned int n, const double eta ){
  double err = 0;

  while( err <= eta && n-- ){
    err += fabs( series[n] );
  }

  return n++;
}

static DDG::Complex
s11( const double t ){
  // how many entries to use?
  static unsigned int ns11r = 0, ns11i = 0;

  if( !ns11r ){
    ns11r = inits( s11r, sizeof s11r / sizeof *s11r, mach[2] / 10 );
    ns11i = inits( s11i, sizeof s11i / sizeof *s11i, mach[2] / 10 );
  }
  return DDG::Complex( csevl( t, s11r, ns11r ), csevl( t, s11i, ns11i ) );
}

static DDG::Complex
s12( const double t ){
  // how many entries to use?
  static unsigned int ns12r = 0, ns12i = 0;

  if( !ns12r ){
    ns12r = inits( s12r, sizeof s12r / sizeof *s12r, mach[2] / 10 );
    ns12i = inits( s12i, sizeof s12i / sizeof *s12i, mach[2] / 10 );
  }
  return DDG::Complex( csevl( t, s12r, ns12r ), csevl( t, s12i, ns12i ) );
}

static DDG::Complex
m12( const double t ){
  // how many entries to use?
  static unsigned int nm12r = 0, nm12i = 0;

  if( !nm12r ){
    nm12r = inits( m12r, sizeof m12r / sizeof *m12r, mach[2] / 10 );
    nm12i = inits( m12i, sizeof m12i / sizeof *m12i, mach[2] / 10 );
  }
  return DDG::Complex( csevl( t, m12r, nm12r ), csevl( t, m12i, nm12i ) );
}

DDG::Complex
DirichletIJDirect( const double s, const double gii, const double gij, const double gjj ){
  const double s2(s*s);
  const DDG::Complex is(0,s), is3(0,s*s2);
  const double s4(s2*s2);
  return
    (((3*gii + 4*gij + 3*gjj) + is*(gii + gij + gjj) - is3*gij/6
      + DDG::Complex(cos(s),sin(s))*(-(3*gii + 4*gij + 3*gjj)
				      + is*(2*gii + 3*gij + 2*gjj)
				      + s2*(gii + 2*gij + gjj)/2))/s4
     + (gii - 2*gij + gjj)/24 - is*(gii - 2*gij + gjj)/60);
}

// gii = |e_{ki}|^2, gij = <e_{ki}, e_{kj}>, gjj = |e_{jk}|^2
DDG::Complex
DirichletIJ( const double s, const double gii, const double gij, const double gjj ){
  if( fabs(s) > M_PI ){
    return DirichletIJDirect( s, gii, gij, gjj );
  }else if( s > 0 ){
    // between 0 and pi
    const double t = s*2/M_PI - 1; // normalize
    return ((gii + gjj)*s11(t) + gij*s12(t)).conj();
  }else{
    // between -pi and 0
    const double t = -s*2/M_PI - 1; // normalize
    return ((gii + gjj)*s11(t) + gij*s12(t));
  }
}

// gjj = |e_{ij}|^2, gjk = <e_{ij}, e_{ik}>, gkk = |e_{ki}|^2
double
DirichletII( const double s, const double gjj, const double gjk, const double gkk ){
  return .25*((gjj-2*gjk+gkk)+s*s*(gjj+gjk+gkk)/90);
}

double
MassII( void ){ return 1/6.; }

DDG::Complex
MassIJDirect( const double s ){
  const double s2 = s*s;
  const double s4 = s2*s2;
  const DDG::Complex is(0,s), is3(0,s*s2);
  return (6*DDG::Complex(cos(s), sin(s)) - 6 - 6*is + 3*s2 + is3)/(3*s4);
}

DDG::Complex
MassIJ( const double s ){
  if( fabs(s) > M_PI ){
    return MassIJDirect( s );
  }else if( s > 0 ){
    const double t = s/M_PI*2-1;
    return m12(t).conj();
  }else{
    const double t = -s/M_PI*2-1;
    return m12(t);
  }
}
